%%Set up Matlab
clear all;
estimate=1;
 
%%Set parameters;
alpha =  0; 
alpha_2 = 0;
kappa = 0; 
lambda_4 = 0;  
lambda_5 = 0;
ln_sigma_v = 0; 
theta =[alpha alpha_2 lambda_4 lambda_5  kappa]'; 
y = [theta' ln_sigma_v]; 
y = y';
x = ["Parameter" "mu_v" "alpha" "S_4" "S_5" "kappa" "sigma_v"];
x = x';
v=0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Estimate the model %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if estimate==1;
    
%%Load data;
  filename = 'dynamic_data_2022.csv';
  R=importdata(filename);
  R=R.data;
      
%%Initial guess;
theta_0 = y;
   
%%Optimization options;
options_fmin = optimoptions('fminunc', 'Display', 'iter', 'Algorithm', 'quasi-newton', ...
   'MaxFunEvals', 10000, 'MaxIter', 500, 'SpecifyObjectiveGradient', false, ...
   'OptimalityTolerance', 1e-6);

%%%Maximize likelihood   
tic
[theta_hat Qhat exitflag output G H]=fminunc(...
    @(x)likelihood_objective(x,R),theta_0,options_fmin);
toc

%%%Display estimates and standard errors;
%Compute standard errors;
SE=sqrt(diag(inv(H)));

C_1=R(:,1);
C_2=R(:,2);
C_3=R(:,3);
C_4=R(:,4);
C_5=R(:,5);
X_1=R(:,6);
D_1=R(:,7);
N=length(C_1);

%Verify momemnts from insheeted empirical data
Pr5    = nanmean(C_5(X_1==0))
Pr5_S4 = nanmean(C_5(X_1==2))
Pr5_S5 = nanmean(C_5(X_1==3))

Pr4    = nanmean(C_4(X_1==0))
Pr4_S4 = nanmean(C_4(X_1==2))
Pr4_S5 = nanmean(C_4(X_1==3))

Pr3    = nanmean(C_3(X_1==0))
Pr3_S4 = nanmean(C_3(X_1==2))
Pr3_S5 = nanmean(C_3(X_1==3))

Pr2    = nanmean(C_2(X_1==0))
Pr2_S4 = nanmean(C_2(X_1==2))
Pr2_S5 = nanmean(C_2(X_1==3))

Pr1    = nanmean(C_1(X_1==0))
Pr1_S4 = nanmean(C_1(X_1==2))
Pr1_S5 = nanmean(C_1(X_1==3))

%Display results;
[theta_hat SE];

theta_hat = theta_hat';
theta_hat = ["Estimate" theta_hat];
theta_hat = theta_hat';

theta_0 = theta_0';
theta_0 = ["Starting value" theta_0];
theta_0 = theta_0';

SE = SE';
SE = ["SE" SE];
SE = SE';

m = [x theta_0 theta_hat SE] 
writematrix(m,'structuralmodel.csv') 

[Pr5  Pr4  Pr3  Pr2  Pr1]
 
Control = [Pr5  Pr4  Pr3  Pr2 Pr1]'
Signal5 = [Pr5_S5  Pr4_S5  Pr3_S5  Pr2_S5 Pr1_S5]'
Signal4 = [Pr5_S4  Pr4_S4  Pr3_S4  Pr2_S4 Pr1_S4]'

Sim = [Control(:), Signal5(:), Signal4(:)];

end;